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EIGENSOLUTION  ANALYSIS  OF  THE  DISCONTINUOUS  GALERKIN  METHOD 
WITH  NON-UNIFORM  GRIDS,  PART  I:  ONE  SPACE  DIMENSION 

FANG  Q.  HU*  AND  HAROLD  L.  ATKINS+ 


Abstract.  We  present  a  detailed  study  of  spatially  propagating  waves  in  a  discontinuous  Galerkin 
scheme  applied  to  a  system  of  linear  hyperbolic  equations.  We  start  with  an  eigensolution  analysis  of  the 
semi-discrete  system  in  one  space  dimension  with  uniform  grids.  It  is  found  that,  for  any  given  order  of 
the  basis  functions,  there  are  at  most  two  spatially  propagating  numerical  wave  modes  for  each  physical 
wave  of  the  Partial  Differential  Equations  (PDE).  One  of  the  modes  can  accurately  represent  the  physical 
wave  of  the  PDE  and  the  other  is  spurious.  The  directions  of  propagation  of  these  two  numerical  modes 
are  opposite,  and,  in  most  practical  cases,  the  spurious  mode  has  a  large  damping  rate.  Furthermore, 
when  an  exact  characteristics  split  flux  formula  is  used,  the  spurious  mode  becomes  non-existent.  For  the 
physically  accurate  mode,  it  is  shown  analytically  that  the  numerical  dispersion  relation  is  accurate  to  order 
2p  +  2  where  p  is  the  highest  order  of  the  basis  polynomials.  The  results  of  eigensolution  analysis  are  then 
utilized  to  study  the  effects  of  a  grid  discontinuity,  caused  by  an  abrupt  change  in  grid  size,  on  the  numerical 
solutions  at  either  side  of  the  interface.  It  is  shown  that,  due  to  “mode  decoupling”,  numerical  reflections  at 
grid  discontinuity,  when  they  occur,  are  always  in  the  form  of  the  spurious  non-physical  mode.  Closed  form 
numerical  reflection  and  transmission  coefficients  are  given  and  analyzed.  Numerical  examples  that  illustrate 
the  analytical  findings  of  the  paper  are  also  presented. 

Key  words,  finite  element  methods,  unstructured  grids,  wave  propagation,  acoustics 

Subject  classification.  Numerical  Analysis 

1.  Introduction.  Discontinuous  Galerkin  Method  (DGM)  is  a  finite  element  method  that  allows  a 
discontinuity  of  the  numerical  solution  at  element  interfaces.  It  has  been  developed  very  rapidly  in  the  past 
few  years  and  has  been  applied  to  many  fields  of  practical  importance,  such  as  computational  fluid  dynamics, 
aeroacoustics,  and  electromagnetics  (e.g.,  [2,3,7,22]).  A  recent  review  of  DGM  can  be  found  in  [5]  with  an 
extensive  list  of  references. 

It  is  well  known  that  for  a  discontinuous  Galerkin  scheme  employing  basis  polynomials  up  to  order  p,  the 
rate  of  convergence  is  in  general  and  in  some  special  cases  ([13,12,17,19]),  where  /i  is  a  measure 

for  the  size  of  elements.  Occurrences  of  super-convergence  in  DGM  has  been  reported  in  the  literature  and 
some  are  reviewed  in  [5]  .  For  examples,  Biswas,  Devine  and  Flaherty  [4]  and  Adjerid,  Aiffa  and  Flaherty  [1] 
showed  super-convergence  on  Gauss-Radau  points.  Lowrie,  Roe  and  van  Leer  reported  numerical  results  of 
order  2p+l  convergence  in  [14].  Most  recently,  Cockburn,  Luskin,  Shu  and  Siili  [6]  showed  the  possibility 
of  obtaining  2p+l  convergence  by  a  suitable  post-processing  of  the  numerical  solution. 

In  contrast  to  the  many  studies  on  convergence  rates,  there  have  been  relatively  fewer  works  on  the 
wave  propagation  properties  of  DGM.  In  [12],  Johnson  and  Pitkaranta  included  a  Fourier  analysis  of  DGM 
for  the  case  of  p  =  1  and  showed  that  the  eigenvalue  of  the  “amplification  matrix”  (E(0))  is  accurate  to 
order  4  (local  error).  In  [15],  Lowrie  performed  a  Fourier  analysis  of  a  space-time  discontinuous  Galerkin 

‘Department  of  Mathematics  and  Statistics,  Old  Dominion  University,  Norfolk,  VA  23529.  Supported  by  the  National 
Aeronautics  and  Space  Administration  under  NASA  Contract  No.  NASl-97046  while  the  first  author  was  in  residence  at 
ICASE,  NASA  Langley  Research  Center,  Hampton,  VA  23681,  and  by  NASA  Langley  Grant  NAG-1-01044. 
tNASA  Langley  Research  Center,  Hampton,  VA  23681. 
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scheme,  up  to  p  =  3,  for  a  one-dimensional  scalar  advection  equation  and  showed  that  the  eigenvalue  is 
accurate  to  order  2p  +  2  (locally)  which  results  in  a  global  order  2p  +  1  decay  of  the  evolution  component 
of  the  numerical  error.  In  [9],  Hu,  Hussaini  and  Rasetarinera  studied  numerical  dissipation  and  dispersion 
errors  of  DGM  for  one-  and  two-dimensional  wave  equations.  They  also  analyzed  anisotropic  errors  of  wave 
propagation  in  triangular  and  quadrilateral  elements.  In  a  recent  work  by  Rasetarinera,  Hussaini  and  Hu 
[18],  it  was  further  demonstrated  numerically  that  dissipation  errors  of  DGM  decay  at  order  2p  +  2  (locally) 
when  the  exact  characteristics  splitting  flux  formula  is  used.  Another  study  of  Fourier  analysis  was  carried 
out  in  [20]  by  Sherwin  which  gave  exact  expressions  of  the  numerical  frequency  analytically  up  to  p  =  3 
and  numerically  for  p  =  10.  It  is  also  interesting  to  note  here  some  related  works  in  continuous  Galerkin 
methods.  For  instance,  a  case  of  super-convergence  in  phase  error  in  continuous  Galerkin  methods  has  been 
shown,  numerically,  by  Thompson  and  Pinsky[21]  and,  theoretically,  by  Ihlenburg  and  Babuska[ll].  It  was 
found  that  when  basis  polynomials  of  order  p  are  used  the  phase  error  converges  (locally)  at  [11].  We 

point  out  that  this  is  one  order  less  than  that  for  DGM  as  we  will  show  in  this  paper. 

The  present  work  has  been  motivated  primarily  by  the  need  to  understand  wave  propagation  in  DGM 
with  non-uniform  elements  (grids).  As  a  first  step  toward  such  a  goal,  we  study  wave  propagation  through 
an  interface  with  an  abrupt  change  in  grid  size  in  one  space  dimension.  We  will  first  carry  out  an  analysis  on 
spatially  propagating  waves,  referred  to  as  the  eigensolutions,  of  the  semi-discrete  system  in  uniform  grids. 
Then  the  results  of  such  an  analysis  will  be  applied  to  study  wave  reflection  and  transmission  by  expressing 
the  numerical  solution  on  either  side  of  the  interface  in  eigensolutions.  The  reflection  and  transmission 
coefficients  are  then  found  by  deriving  proper  coupling  conditions  at  the  interface.  As  we  will  see,  numerical 
reflection  at  a  grid  discontinuity  is  dependent  on  the  flux  formula  employed  in  the  implementation  of  DGM. 
Two  commonly  used  flux  schemes  are  considered  in  this  paper,  namely,  the  characteristics-based  flux  and 
Lax-Friedrich  flux  formulas.  These  two  schemes  will  be  analyzed  in  a  unified  way  by  introducing  an  upwind 
factor. 

We  point  out  that  a  major  difference  between  the  present  and  previous  works  in  wave  analysis  for 
DGM  is  that  in  the  present  work  we  study  spatial  waves  where  the  temporal  frequency  is  specified  and 
the  corresponding  wavenumber  is  sought  as  eigenvalues,  while  in  the  previous  studies  the  wavenumber  was 
specified  and  the  frequency  was  found  as  eigenvalues.  The  present  approach  is  necessary  because,  with 
an  introduction  of  grid  discontinuity,  numerical  wave  number  is  not  constant  across  the  interface  of  a  grid 
change.  The  use  of  spatial  waves  also  turn  out  to  be  advantageous  in  that  the  eigenvalue  problem  is  greatly 
simplified  and  reduced.  As  a  result,  the  numerical  dispersion  relation  is  governed  by  a  quadratic  equation 
that  can  be  solved  analytically  for  any  order  of  the  basis  polynomials.  Specifically,  in  a  uniform  grid,  it 
is  found  that  there  are  at  most  two  spatially  propagating  numerical  wave  modes  for  each  physical  wave  of 
the  PDF.  One  of  the  numerical  wave  modes  can  accurately  approximate  the  physical  wave  and  the  other 
is  a  highly  irregular  spurious  mode.  They  will  be  referred  to  as  the  physical  and  spurious  numerical  waves 
respectively  in  this  paper.  For  the  physically  accurate  mode,  it  will  be  shown  that  the  numerical  dispersion 
relation  is  accurate  to  (A:/i)^^+^  locally  where  k  is  the  wave  number,  which  confirms  those  previous  works 
mentioned  earlier[12][14][18].  In  fact,  we  will  show  that  dispersion  error  is  of  order  2p-|-3  while  the  dissipation 
error  is  of  order  2p  +  2.  For  the  spurious  mode,  it  is  found  that  it  propagates  in  the  opposite  direction  of  the 
physical  mode  and  becomes  non-existent  when  the  exact  characteristics-based  flux  formula  is  used.  Following 
the  analysis  of  waves  in  uniform  grids,  the  effect  of  a  grid  change  on  either  side  of  the  interface  is  studied.  It 
is  found  that  waves  associated  with  different  physical  eigenvectors  are  decoupled  and  numerical  reflections 
are  always  in  the  form  of  the  spurious  numerical  wave  and  are  highly  damped. 
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The  rest  of  the  paper  is  organized  as  follows.  In  section  2,  we  describe  the  discretization  process  and  the 
associated  flux  formulas.  In  section  3,  the  eigenvalue  problem  for  spatially  propagating  waves  in  a  uniform 
grid  is  formulated.  In  section  4,  numerical  dispersion  relation  and  its  accuracy  are  analyzed  and  discussed. 
Wave  propagation  through  a  grid  discontinuity  is  studied  in  section  5.  And  numerical  examples  are  presented 
in  section  6.  Section  7  has  our  conclusions. 

2.  Formulations  of  discretization  and  numerical  flux.  Consider  the  discontinuous  Galerkin  method 
for  a  system  of  hyperbolic  equations  in  one-dimensional  space: 


(2.1) 


3u  3f(u) 
dt  dx 


where  u  is  a  vector  of  dimension  N  and  f  is  the  flux  vector.  We  will  only  consider  linear  cases  in  our  analysis 
and  assume  that 


(2.2)  f(u)  =  An 

where  A  is  a  constant  N  x  N  matrix.  We  assume  that  A  has  N  real  eigenvalues,  denoted  by  aj  for 
j  =  1,2,  ...,N  and  the  eigenvectors  of  A,  denoted  by  ej,  form  a  complete  basis  in  A-dimensional  space. 
Throughout  this  paper,  unless  specified  otherwise,  lower  case  bold  face  letters  will  stand  for  column  vectors 
and  upper  case  bold  face  letters  stand  for  matrices. 

In  a  discretization  of  (2.1)  using  the  discontinuous  Galerkin  method,  the  spatial  domain  is  partitioned 
into  elements,  =[xn-i ,  Xn],  where  n  is  the  element  index.  In  each  element,  the  numerical  solution,  denoted 
by  vL^{x,t),  is  expressed  as 


(2.3)  Uft(a;,t)  =  ^c'^{t)p'l{x), 

where  {p^{x),<i  =  0, 1,  ...,p}  is  the  set  of  basis  polynomials  for  element  i?„.  Here  p,  without  superscript  or 
subscript,  denotes  the  highest  order  of  polynomials  in  the  chosen  basis  and  c"(t)  is  the  expansion  coefficient. 
In  a  weak  formulation  for  (2.1),  we  require  that 

,2.4)  +  =  0 

for  a'  =  0, 1,  ...,p.  By  a  use  of  integration  by  parts,  the  above  is  re-written  as  follows. 


(2.5) 


dvil 


■p^,(x)dx+  [f^-p^,(x)]: 


ox 


At  any  interface  between  two  elements,  i.e.,  the  end  points  Xn-i  and  Xn,  the  flux  vector  f-^  is  not 
uniquely  determined  and  a  flux  formula  has  to  be  supplied  to  complete  the  discretization  process.  Various 
kinds  of  flux  formulas  have  been  proposed  and  used  in  the  literature.  In  this  paper,  we  will  consider  two 
commonly  used  flux  formulas.  They  are  specified  below  and  will  be  referred  to  as  the  characteristics-based 
flux  formula  and  Lax-Friedrich  flux  formula  respectively. 

The  characteristics-based  flux  formula  is  of  the  form 


3 


(2.6) 


6»  >  0 


f^{uL,un)  =  ^[f(uL)  +f(uij)  -6»|A|(uij  -u^)], 

where  and  are  the  values  of  u  at  the  interface  calculated  using  expansion  coefficients  of  the  elements  at 
the  left  and  right  of  that  interface  respectively.  (In  DGM  and  are  not  required  to  be  the  same.)  Here 
0  is  a  scalar  parameter.  The  value  of  0  is  usually  unity  in  practice  which  makes  (2.6)  an  exact  characteristics 
splitting  (the  exact  Roe  solver).  On  the  other  hand,  (2.6)  will  result  in  a  symmetric  averaged  scheme  when 
6  =  0.  Here,  we  will  keep  0  as  a  parameter  so  that  our  analysis  can  be  useful  for  a  wide  range  of  cases.  For 
convenience,  (2.6)  will  be  written  as 


(2.7)  f^(uL,  Uij)  =  Alul  +  Arur 

where 


(2.8)  Ai  =  i[A+0|A|],  A^  =  i[A-0|A|]. 

The  Lax-Friedrich  flux  formula  is  of  the  form 

(2.9)  f^(uL,Uij)  =  i[f(uL) +f(uij) -6»|a|ma;j(uij -ul)]  6>0 

where  \a\max  is  the  maximum  (absolute  value)  of  the  eigenvalues  of  A.  This  can  again  be  written  in  the 
form  of  (2.7)  with 

(2.10)  Ai  =  -[A  +  6*|a|maa;I],  =  -  [A  -  6*|a|maa:I]  • 

Using  expression  (2.7)  for  both  cases,  the  semi-discrete  equation  (2.5)  can  now  be  written  as 

J  +  [A.Lul{x„,t)  +  ARul+^{x„,t)]  Pi,{x„) 


(2.11)  -  [ALu”“^(a;„_i,t)  +  Aiju;((a;„_i,t)]  p”,(a;„_i)  -  y  Au)(-^da;  =  0 

for  P  =  0, 1, 

Together  with  (2.3),  equation  (2.11)  yields  a  system  of  time  evolution  equations  for  the  expansion 
coefficients  for  each  element.  This  system  is  usually  solved  by  some  time  integration  scheme  such  as  the 
Runge-Kutta  schemes  ([2],  [7]). 
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3.  Spatially  propagating  waves  in  nniform  grids. 

3.1.  Use  of  local  variables.  Introducing  a  local  coordinate  ^  for  each  element,  we  let 


(3.1)  ^  - (x-Xn)  where  Ax„  =  x„  -  x„-i  and  Xn  =  — — - -■ 

Axn  2 

In  addition,  the  basis  functions  will  be  taken  to  be  the  same  for  all  the  elements  when  expressed  in  the  local 

coordinate  i.e.,  we  assume 

Peix)  =  PtiO, 

where  {Pi{Ot^  =  Oj  ■■■tP}  is  a  chosen  set  of  basis  functions,  such  as  the  Legendre  polynomials,  or  the  set 
of  {1,^,^^,  The  results  of  our  analysis  are  independent  of  the  specific  choice  on  the  basis  functions. 

We  look  for  wave-like  solutions  supported  by  (2.11).  By  assuming  a  periodicity  in  time  with  a  frequency 
Lo,  we  let 


(3.2)  u”(^,  t)  =  e--*u”(a  where  uKO  =  Y.  and  *  = 

e^o 

The  expansion  coefficients  c"  are  now  independent  of  t.  Substituting  the  above  into  (2.11),  we  get 

111}  /\  'IT  A ^ 

- ^  J  1  +  [Al<(1)  +  A^u”+i(-l)]  PHI) 

(3.3)  -  [AiuHi(l)  +  A^uH-1)]  H'(-l)  -  AuliO^d^  =  0 
for  £'  =  0, 1,  ...,p. 

3.2.  Uniform  grid  and  the  eigenvalne  problem.  We  now  consider  the  case  where  elements  are 
uniform  in  length,  i.e.,  Axn  =  h.  After  substituting  (2.3)  into  (3.3),  we  look  for  solutions  with  expansion 
coefficients  of  the  form 


(3.4)  H  =  A”c^ 

where  A  is  an  undetermined  complex  number  and  ci  is  a  vector  independent  of  the  element  index  n.  It  is 
easy  to  see  that,  if  we  express  A  as 

(3.5)  A  =  e**'*'', 

then,  kh  can  be  interpreted  as  the  wave  number  of  the  numerical  solution.  Here  kh  will  be  referred  to  as  the 
numerical  wave  number. 

For  convenience  of  discussion,  we  define  a  column  vector  that  contains  all  the  expansion  coefficients 


(3.6) 


Co 

Cl 
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and  matrices 


(3.7)  Q  =  {g'f4>  where  qi'i  =  J  PiiOPfiOd^, 

(3.8)  Q'  =  {g^,  where  =  J  ^  PiiO^^d^, 

(3.9)  =  {bi>i},  where  hti  =  Pi<{a)Pi{h) 
where  £  =  0, 1,  ...,p.  Then,  equation  (3.3)  can  be  re-written  compactly  as 

-^(Q  ®  I)*  +  (B(1,1)  ®  Al)x  -  (B(_1  _1)  ®  Aij)x  -  (Q'  ®  A)x 


(3.10)  +A(B(i  _i)  g)  Aij)x  -  ^(B(_i_i)  g)  Al)x  =  0 

where  g)  denotes  the  Kronecker  product  (The  definition  and  relevant  properties  of  g)  can  be  found  in  appendix 
Al).  For  a  given  frequency  to,  equation  (3.10)  forms  an  eigenvalue  problem  with  A  being  the  eigenvalue  and 
X  the  eigenvector. 

Next,  we  show  that  (3.10)  can  be  equivalently  separated  into  N  independent  eigenvalue  problems,  where 
each  sub-problem  corresponds  to  one  of  the  physical  wave  modes  of  the  PDE.  Here  each  pair  of  the  eigenvalue 
and  eigenvector  {aj,ej}  of  the  PDE  will  be  referred  to  as  a  wave  mode  of  the  PDE  (2.1)  and  aj  is  the  wave 
speed  of  that  mode.  Since  we  are  interested  in  spatially  propagating  waves,  we  assume  aj  ^  0. 

We  first  express  x  given  in  (3.6)  in  terms  of  eigenvectors  of  the  PDE  as  follows. 


Co 

1 

M 

it  2 

o 

_ 1 

yoj 

(3.11) 

X  = 

Cl 

= 

yii®i 

II 

yi3 

3  =  i 

Cp 

yp3^3 

ypj 

N 


where  yj  is  a  column  vector  of  dimension  p  -I-  1.  By  substituting  (3.11)  into  (3.10)  and  using  a  property  of 
the  Kronecker  product  (equation  (7.2)  of  Appendix  Al),  we  get 


®  +  (B(i.i)yi)  ®  A^ej  -  (B(_i  _i)yj)  ®  A^ej  -  (Q'yj)  ®  Ae^ 

i=i  ^ 


(3.12) 


-|-A(B(i  _i)yj)  g)  Anej 


T(B(_i,i)yj)  ®  A^e, 


=  0. 


Furthermore,  for  Al  and  A^  given  in  (2.8)  and  (2.10),  we  have 


(3.13) 


Aie,  =  i±^a,e 


3  ^3  ■> 


— 


i-M 


^3 '^3, 
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where 


(3.14)  /J,  =  M 

aj 

for  the  characteristics-based  flux  deflned  in  (2.8)  and 


(3.15)  I3j  = 

Uj 

for  the  Lax-Friedrich  flux  deflned  in  (2.10).  Thus  the  two  types  of  flux  formulas  can  be  treated  in  a  unified 
way  by  using  (3.13).  For  convenience  of  discussion,  we  define 


(3.16)  7j=dpj. 

and  7j  will  be  referred  to  as  the  upwind  factor  of  the  scheme  for  the  jth  wave  mode  of  the  PDF.  We  note 
that,  for  both  cases  given  in  (3.14)  and  (3.15),  \^j\  =  1  leads  to  the  exact  characteristics  splitting.  For 
|7j|  >  1,  the  eigenvalues  of  and  are  all  positive  and  negative,  respectively.  Note  also  that  for  the 
Lax-Friedrich  flux  formula  applied  to  a  system  of  equations  in  which  the  \aj\  varies  widely,  |7j|  will  be  large 
for  the  slowest  wave  modes;  however,  the  wavenumber  Loh/oj  will  also  be  proportionally  large  for  any  given 

CJ. 

Equation  (3.12)  can  now  be  expressed  as 

N 

El  iioh  1  +  7?  1  “  7?  r^i 

I - «iB(i,i)yj - -  OjQ  yj 

(3.17)  _i)yj  -  j  ®  ej  =  0. 

Due  to  linear  independency  of  e^’s,  it  is  easy  to  see  that  (3.17)  yields  N  independent  sub-eigenvalue 
problems. 


(3.18) 


-  — Q  +  (1  +  7i)B(i,i)  -  (1  -  7i)B(-i,-i)  -  2Q'  +  A(1  -  7i)B(i,_i)  -Ul  +  7i)B(_i,i)l  y^ 

Clj  A  . 


=  0 


for  j  =  1,  2, ...,  A,  in  which  yj  is  the  eigenvector  and  A  is  the  eigenvalue. 

We  observe  that,  by  solving  the  eigenvalue  problem  posed  in  (3.18),  we  will  obtain  A  as  a  function  of 
the  non-dimensional  frequency  ^  (or  wavenumber)  and  the  upwind  factor  7^,  i.e.. 


(3.19)  A  =  F(— ,7i). 

Oj 

Since  A  is  directly  related  to  the  numerical  wave  number  kh  by  (3.5),  equation  (3.19)  is  the  numerical 
dispersion  relation  of  the  scheme.  It  is  an  intrinsic  property  of  the  discretization. 

In  addition,  the  non-trivial  solution  of  (3.18)  forms  the  eigenfunction  of  the  numerical  mode.  Specifically, 
let  the  eigenvectors  of  (3.18)  be  denoted  by  y  =  {n^},  then  the  eigenfunctions  will  be  of  the  form 
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(3.20)  uliU)  =  , 

where 


(3.21)  <(a=e“^'*V(^;— ,7i),  and  /(^;  ^^.)  =  ^  ^;,P,(^) 

Gi  'i  G  'i 

3  3 

For  convenience,  the  eigenfunctions  will  be  normalized  such  that 

(3.22)  = 

ttj 

4.  Numerical  dispersion  relation. 

4.1.  Determinant  of  equation  (3.18).  Wave  propagation  properties  of  the  numerical  scheme  are 
encoded  in  the  numerical  dispersion  relation  (3.19).  For  convenience  of  discussion,  the  subscript  j  will  be 
dropped  in  this  section.  We  define 

K  =  —  and  =  k^h, 

a 

where  K  is  the  non-dimensional  exact  wave  number  of  the  PDF  and  is  the  non-dimensional  numerical 
wave  number  as  given  in  (3.5).  is  related  to  A  of  (3.18)  by 

(4.1)  A  =  e*'^7 

By  letting  the  determinant  of  the  coefficient  matrix  for  yj  in  (3.18)  be  zero,  we  get  an  algebraic  equation 
for  A  for  any  given  value  of  K.  We  have  computed  the  determinant  of  (3.18)  symbolically  using  the  computer 
algebra  system  MAPLE  [16].  It  is  found  that  the  determinant,  after  some  normalization,  can  be  written  in 
the  following  form: 


(4.2)  (1  -  -f)[G{iK)X  -  H{iK)]  +  (-l)^’+i(l  +  =  0 

A 

where  G{x)  and  H{x)  are  polynomials  of  degree  p  and  p+  1,  respectively,  with  real  coefficients.  The  exact 
expressions  for  G  and  H  are  given  in  appendix  A2. 

Before  presenting  the  numerical  and  analytical  results  of  (4.2),  we  make  some  general  remarks: 

1.  We  observe  that,  due  to  the  fact  that  the  rank  of  B  matrices  in  (3.18)  is  unity  (see  equation  (3.9)), 
equation  (4.2)  is  quadratic  in  A.  Consequently,  there  will  be  at  most  two  distinct  solutions  for  A  in  (4.2). 
Furthermore,  when  7  =  ±1  (exact  characteristics  splitting  of  the  flux),  equation  (4.2)  becomes  linear  in  A 
and  there  will  be  only  one  solution  for  A. 

2.  In  previous  works  in  the  literature,  the  spatial  wave  number  K^,  thus  A  in  (4.2),  is  specified  and  K  is 
to  be  solved  from  (4.2).  This,  of  course,  will  result  in  solving  a  polynomial  of  degree  p+1  which  is  difficult 
to  do  analytically  for  p  >  3. 

3.  If  A  =  Ao  is  a  solution  to  (4.2)  for  7  =  70,  then  A  =  ^  is  a  solution  for  7  =  —70.  (This  can  be  shown 
simply  by  taking  a  complex  conjugate  of  (4.2).)  Thus,  it  is  sufficient  to  consider  theoretically  only  the  cases 
with  7  >  0  in  (4.2),  i.e.,  right-going  waves,  for  dissipation  and  dispersion  errors. 
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4.2.  Numerical  results  of  the  eigensolutions.  We  first  present  numerical  results  of  (4.2).  Its 
analytical  properties  will  be  presented  in  section  4.3.  As  we  have  seen  in  previous  discussions,  when  |7|  ^  1, 
there  are  two  roots  for  A.  Each  root  represents  a  numerical  wave  mode  whose  wave  number  is  found  by 
(4.1)  and  whose  mode  shape  (eigenfunctions)  found  by  (3.21).  As  we  will  see,  one  of  the  numerical  modes 
can  faithfully  represent  the  physical  wave  and  the  other  mode  is  spurious  or  non-physical.  These  two  modes 
behave  very  differently  and  it  is  easy  to  distinguish  the  physical  mode  from  the  non-physical  mode.  When 
I7I  =  1  (exact  characteristics  flux),  of  course,  there  will  be  only  one  root  and  the  spurious  mode  will  not  be 
present. 

We  will  use  a  case  with  7  =  0.5  as  an  example  to  demonstrate  numerical  results.  For  a  given  value  of 
exact  wave  number  K,  we  solve  equation  (4.2)  and  obtain  two  values  of  A  which  are  then  converted  into 
numerical  wave  numbers  according  to  (4.1).  For  the  physical  mode,  the  numerical  wave  number  as 
a  function  of  K  is  plotted  in  Figure  4.1(a)  (real  part)  and  4.1(b)  (imaginary  part),  for  cases  p  =  1,2,  3, 4. 
The  diagonal  line  in  Figure  4.1(a)  is  the  exact  relation,  i.e.,  when  the  numerical  wave  number  is  equal  to 
the  actual  physical  wave  number.  It  is  seen  that,  for  a  given  value  of  p,  the  real  part  of  the  numerical 
wave  number,  Re{Kh),  follows  the  exact  line  closely  for  a  range  of  K  values.  This  range  will  be  termed 
resolved  wave  number  spaee.  Clearly,  the  higher  the  order  of  the  basis  functions,  the  larger  the  resolved 
space.  Figure  4.1(b)  shows  the  imaginary  part  of  the  numerical  wave  number,  Im{Kh).  Since  the  wave  is 
right-traveling  for  the  present  case  (7  >  0) ,  the  positive  imaginary  part  represents  numerical  damping  as  the 
wave  propagates  in  space.  We  note  that  the  damping  is  not  significant  for  wave  numbers  within  the  resolved 
wave  number  space  in  each  scheme.  The  exact  boundary  of  resolved  range  is,  of  course,  somewhat  arbitrary 
and  depends  on  the  accuracy  criteria  imposed.  This  issue  will  be  closely  examined  in  section  5.2.  In  general, 
the  dissipation  error  places  a  higher  requirement  on  the  resolution  of  the  scheme  than  the  dispersion  error 
in  DGM. 

For  the  spurious  mode,  the  relation  of  Kh  v.s.  K  is  plotted  in  Figure  4.2.  For  the  real  part  of  Kh  shown 
in  Figure  4.2(a),  the  curve  starts  at  0  for  p  odd  and  starts  at  tt  for  p  even.  The  group  velocity  of  these  waves 
(slope  of  Re{Kh)  v.s.  K)  is  negative,  indicating  that  the  spurious  waves  are  left-traveling,  in  the  opposite 
direction  of  the  actual  physical  wave.  The  imaginary  part  of  Kh  is  also  negative,  indicating  again  that  the 
wave  is  left-traveling  and  damped.  The  damping  rates  for  the  spurious  modes  in  Figure  4.2(b)  are  quite 
large  for  the  cases  shown.  This  means  that  the  spurious  mode  is  expected  to  be  damped  very  rapidly  in 
computation. 

The  corresponding  eigenfunctions  of  the  physical  and  spurious  modes  are  plotted  in  Figures  4.3  and  4.4 
respectively.  The  eigenfunctions  are  constructed  according  to  (3.21)  using  eigenvectors  from  (3.18)  as  the 
expansion  coefficients.  Plotted  are  eigenfunctions  over  a  span  of  30  elements,  with  the  first  element  being 
[-1,1]  as  indicated  by  dark  lines  in  the  plots.  As  shown  in  Figure  4.3,  the  physical  mode  travels  to  the  right 
and  the  amount  of  damping  is  quite  visible  for  p  =  1  and  2  with  the  chosen  value  oi  K  =  2.  The  damping 
error  reduces  significantly  as  order  increases. 

In  Figure  4.4,  we  see  that  the  spurious  non-physical  mode  is  damped  very  rapidly  for  all  cases  shown, 
which  is  consistent  with  our  observation  in  Figure  4.2. 

As  will  be  shown  later  (equation  (40)),  the  damping  factor  of  the  spurious  mode  is  related  to  the  value 
of  7  as  plotted  in  Figure  4.5.  Thus  the  spurious  wave  modes  become  highly  damped  when  7  is  close 

to  unity  and  much  less  damped  when  7  is  close  to  zero  or  much  greater  than  unity.  In  practice,  small  7  is 
avoided  by  choosing  0  ~  1;  large  7  occurs  for  slow  wave  modes  where  \a\maxlaj  is  small. 
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K(=wh/a) 


K(=a;h/a) 


Fig.  4.1.  The  physical  mode,  numerical  wave  number  ky,h  v.s.  actual  wave  number  (normalized  frequency)  K  = 

7  =  0.5. 

4.3.  Super-accuracy  of  the  numerical  wave  number.  The  numerical  wave  number  Kf,  of  the 
physical  mode  should  be  a  close  approximation  of  the  actual  wave  number  K,  especially  in  the  long  wavelength 
limit,  i.e.,  when  K  is  small.  Here,  we  give  an  estimation  on  the  order  of  convergence  in  wave  number  space 
and  show  that  Kh  is  accurate  to  the  actual  wave  number  K  to  order  2p  +  2,  which  is  twice  the  order  of 
accuracy  of  the  basis  functions. 

Assuming  |7|  ^  1,  we  can  re-write  (4.2)  as  a  quadratic  equation  for  A  as  follows. 


'  ■*'  1-7  G(iA-)  J '  I--,  GiiK) 

By  examining  the  computed  H{x)  and  G{x)  functions  (appendix  A2),  we  found  that  the  ratio  H{x)/G{x) 
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K=(u;h/a) 


K=(u;h/a) 


Fig.  4.2.  The  spurious  mode,  numerieal  wave  number  kf^h  v.s.  aetual  wave  number  (normalized  frequeney)  K  = 

7  =  0.5. 

is  always  exactly  the  Fade  approximation  of  e*  to  order  2p  +  2.  (This  has  been  calculated  and  verified 
symbolically  for  p  up  to  16  and  is  conjectured  to  be  true  for  all  p.)  That  is,  we  have 

(4.4)  + 

Consequently,  we  can  show  that,  for  K  small,  the  two  roots  of  (4.3)  are 

(4.5)  +  •  •  •  {physical  mode) 
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Fig.  4.3.  Eigen-functions  of  the  physical  mode.  K  =  2.  7  =  0.5.  (a)  p  =  1,  (b)  p  =  2,  (c)  p  =  3,  (d)  p  =  4. 


(4.6)  =  (— 1  +  7  G{  iK)  ^  .  (non  —  physical  mode) 

1  —  7  G(iK) 

where  Ci,  C2  and  L>i,  D2  are  real  coefficients,  and  dots  represent  higher  order  terms  in  (iK).  A  detailed 
derivation  is  given  in  the  appendix  A3.  Here  the  superscripts  (p)  and  (s)  indicate  the  physical  and  spurious 
modes  respectively. 
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7(=e/3) 


Fig.  4.5.  Effects  of  ^  on  the  damping  rate  of  the  spurious  mode. 


(4.7)  =K+  {-l)PiCiK‘^P+‘^  +  {-1)p+'^C!2K‘^p+^  +  ■■■ 

where  Cj  is  also  a  real  coefficient.  Furthermore,  by  considering  the  real  and  imaginary  parts  of  (4.7),  we  can 
get  an  estimation  on  the  convergence  rates  of  the  dispersion  and  dissipation  errors.  Specifically,  we  have 


(4.8)  dispersion  error  :  Re{Kl^^)  -  K  =  {-1)p+^C2K‘^p+^  +  ■  ■  ■ , 


(4.9)  dissipation  error  :  Im(Kj^^)  =  (—1)pCiK‘^p~^‘^  +  ■  ■  ■ . 

That  is,  for  DGM,  the  dominant  error  is  the  dissipation  error  which  reduces  locally  at  order  2p  +  2.  The 
dispersion  error,  on  the  other  hand,  reduces  locally  at  order  2p  +  3.  This  is  confirmed  in  Figure  4.6  where 
the  numerical  dispersion  relations  shown  in  Figure  4.1  are  re-plotted  in  log-log  scale. 

Note  that  when  |7|  =  1,  the  spurious  mode  is  non-existent  and  it  is  straightforward  to  verify  directly 
from  (4.2)  that  (4.7)-(4.9)  are  still  true  for  the  physically  accurate  mode. 

We  also  note  that  for  polynomials  H{x)  and  G{x)  with  given  orders,  (4.4)  is  the  best  possible  order  of 
approximation.  This  suggests  that  (4.7)  is  the  best  asymptotic  numerical  dispersion  relation  possible. 
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order  2p+3 


K(=a;h/a) 


Fig.  4.6.  Local  order  of  convergence  for  the  dispersion  and  dissipation  errors  of  the  physical  mode.  Circle:  numerical 
wave  number  computed  by  equation  (f.2);  solid  line:  theoretical  convergence  rate,  (a)  dispersion  error;  (b)  dissipation  error. 
7  =  0.5. 

5.  Wave  reflection  at  an  interface  of  mesh  discontinnity. 

5.1.  Reflected  and  transmitted  waves.  In  this  section,  we  consider  a  situation  where  the  size  of 
the  element  is  abruptly  changed  from  hi  to  across  the  interface  between  elements  n  =  0  and  n  =  1,  as 
shown  in  Figure  5.1.  We  will  study  the  wave  reflection  and  transmission  at  the  interface.  Specifically,  we  will 
introduce  an  incident  physical  wave,  traveling  from  left  to  right,  and  look  for  the  reflected  and  transmitted 
waves  caused  by  the  grid  discontinuity. 

Using  the  eigenfunction  expression  (3.20),  we  can  express  the  incident  wave  as 
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incident  wave 


Element  index  n=  -3-2-10  1  2  3 


X_3  X_2 


X2  X3 


Fig.  5.1.  A  schematic  of  grid  change. 


(5.1)  Uincident  =  Aoe  , 7io)eio 

«io 

where  Aq  is  the  wave  amplitude.  Here  denotes  the  eigenvector  of  a  right-going  wave  mode  of  the  PDE 

(2.1) ,  and  and  ^^>7io)  1^®  numerical  wave  number  and  eigenfunction  for  that  wave  mode 

found  assuming  a  uniform  mesh  hi.  In  other  words,  (5.1)  satisfies  the  time  harmonic  semi-discrete  equation 
(3.3)  if  Axn  =  hi  is  held  for  all  n.  The  superscript  (p)  in  (5.1)  denotes  that  the  incident  numerical  mode  is 
a  physical  wave  mode.  Likewise,  a  superscript  (s)  will  be  used  to  denote  the  spurious  modes. 

Due  to  the  discontinuity  in  mesh  size,  there  will  be  reflections  at  the  interface.  For  convenience  of 
discussion,  let  and  denote  the  time  independent  solutions  in  the  left  and  right  half-domains  on 

either  side  of  the  interface  respectively.  By  making  a  use  of  (3.21),  we  get 


(p) 


f(p)l 

■'io 


(jjhi 


7io)« 


incident 


+  ^,7io)eio  +  (^;  ^’7i)« 

JJ^JO 


(•) 


Lohi 

^io 


reflected 


(5.2) 


and 


'  n  _ 
right 


(P) 


ip)( 

io 


Loh^ 


>7io)eio  +  5] 


inK 


^  -3^ 
i¥=ia 


+  .  -  w/l2 


7i)« 


transmitted 


(5.3) 

Here,  and  At  are  amplitudes  of  the  reflected  and  transmitted  waves  associated  with  the  wave  and 
Bj  and  Aj  are  those  associated  with  the  other  waves  of  the  PDE.  The  superscripts  +  and  —  in  the  terms 
inside  the  summations  of  (5.2)  and  (5.3)  denote  the  direction  of  propagation  (right-traveling  and  left-traveling 
respectively).  It  will  be  shown  next,  however,  that  all  Bj  and  Aj  are  zero. 


16 


5.2.  Matching  conditions  at  the  interface.  To  derive  matching  conditions  at  the  interface,  we  first 
note  that  iifgff  and  satisfy  equation  (3.3)  for  a  uniform  element  size  hi  and  /i2  respectively.  The 

coupling  of  the  solutions  can  be  found  by  applying  (3.3)  at  the  two  adjacent  elements  near  the  interface  of 
the  grid  discontinuity,  namely,  at  elements  n  =  0  and  n  =  1,  Figure  5.1.  Thus,  from  (3.3),  we  have 
n  =  0: 

f_^  ■  Pi'im  +  PAl) 

(5.4)  -  [Alu-},(1)  +  A^nO  ^,(-1)]  ^,-(-1)  -  =  0 

n  =  1: 

+  l^LKighA)  +  Ai?U?,,w(-l)]  ^^'(1) 

(5.5)  -  [ALulftil)  +  AnAighti-'^)]  Pi'A'^)  -  J  =  0 

These  two  conditions  can  be  simplified  when  we  recognize  the  fact  that  (5.4)  and  (5.5)  will  still  be  true 
when  Aighti~^)  (5.4)  is  replaced  by  1)  and  u[’g^j(l)  in  (5.5)  is  replaced  by  it-^ighti^)  i'O 

reason  stated  at  the  beginning  of  the  section.  Consequently,  the  matching  conditions  (5.4)  and  (5.5)  are 
equivalent  to  the  following  two  equations  that  are  much  more  compact: 


(5.6) 

and 


-^R^lefti  -^R^riahti 


(5.7)  AluIAI)  =  ALU^gAl) 

Now  by  substituting  (5.2)  and  (5.3)  into  (5.6)  and  (5.7),  and  recalling  (3.13),  we  easily  get 

^jo  ^  ^jo  ^ 


,  •sr^  o  iK~  .  V,-/  ,  1  —  7j 

+  2^Bje  (-1;— ,7i)^— ' 


i¥=io 


(5.8)  =  {-I-  —  ,7io)^ 

«i0  ^  ^ 


a^ej 


and 


^  ^  a  e 


2  '^Jo'=Jo  ^  2  '^Jo'=Jo  ^  2 

i¥=ia 


1  +7i 
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(5.9) 


—  2  ^  2-^  2 

i¥=io 

(In  (5.9),  we  have  used  the  fact  that  eigenfunctions  are  normalized  such  that  f\^=i  =  1,  as  in  (3.22).)  Since 
Gj’s  are  linearly  independent,  it  follows  that 


(5.10) 


~  “0,  j  ^  jo- 


This  means  that  no  component  of  wave  modes  other  than  that  of  the  incident  mode  will  be  present  in 
the  reflected  and  transmitted  waves.  This  also  suggests  that  the  reflected  wave  can  only  be  in  the  form  of 
the  spurious  mode,  the  only  opposite  traveling  numerical  wave  for  the  mode.  An  interesting  consequence 
of  this  is  that  when  the  exact  characteristics  splitting  flux  formula  is  used,  there  will  be  no  reflected  wave 
because  the  opposite-traveling  spurious  wave  is  non-existent. 

Further,  by  equaling  the  coefficients  of  in  (5.8)  and  (5.9),  and  assuming  |  1,  we  get  two  coupled 

equations  for  and  At, 


^>7io)  +  A,e*-'-'o/r(-l;  ^>7io)  =  /io (-1;  ^,7io) 


(P) 


r(p)i 


Lohi 


iK 


(•) 


Lohi 


iK 


(P) 


w/lo 


Aq  +  Ar  —  At- 


Solving  the  above,  we  get  the  following  closed  expressions  for  the  reflection  and  transmission  coefficients 


iK 


(P) 


(5.11) 


A,  e  '-V.^°/r(-l;  ^,7io)  -  (-1;  ^,7io) 

^  e<’.o/W(-l;  -  e<’.o/W(-l; 


(p) 


(5.12) 


— 

^0 


4”  (- 1 ;  Sr .  V. )  -  4"’  <-  r  Sr  •  V. ) 


e'*;-,’ 4'l  {- 1;  ^ .  7„ )  -  e<-..  /<f  { - 1;  ^ .  7,. ) 


Thus,  numerical  reflection  and  transmission  coefficients  are  directly  related  to  the  change  in  dispersion 
properties  of  the  scheme  when  grid  change  occurs. 

To  express  the  above  in  a  more  compact  and,  perhaps,  more  insightful  form,  we  note  the  fact  that 
numerical  solutions  in  DGM  have  a  small  discontinuity  (or  gap)  at  the  boundary  of  any  two  elements.  This 
discontinuity,  of  course,  becomes  diminished  with  the  increase  of  the  resolution  of  the  scheme.  Specifically, 
if  we  let  A/i  denote  the  discontinuity  of  the  numerical  solution  at  the  interface  of  elements  n  =  0  and  n  =  1 
had  the  grid  size  been  uniformly  h,  then  we  have 


Aft  =  4j,(-l)  -<(1)  =  e*'^'*/(-l;— ,7)  -  1 

where  u^{^)  is  the  eigenfunction  specified  in  (3.21).  Thus,  in  terms  of  Aft,  the  expressions  for  the  reflection 
and  transmission  coefficients  given  in  (5.11)  and  (5.12)  can  now  be  written  as 
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h2/hi 


Fig.  5.2.  Solid  lines,  boundaries  of  2%  numerieal  refleetion.  Refleeted  wave  is  less  than  2%  of  the  ineident  wave  for 
parameters  above  the  eurves.  7  =  0.5.  Dashed  lines,  aeeuraey  limits  determined  from  the  dispersion  relation  of  a  uniform  grid 
h2-  The  aeeuraey  limit  for  p  =  1  (not  shown)  is  far  above  and  out  of  the  pieture. 


(5.13) 


^0 
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-Ag 
-  Ag 


(5.14) 


^0 


Ag 

Ag 


-Ag 
-  Ag 


in  which  the  superscript  denotes  the  mode  type,  the  physical  (p)  or  spurious  (s)  mode,  and  the  subscript 
denotes  the  mesh  spacing  used  for  calculating  the  gap. 

Equation  (5.13)  implies  that  numerical  reflection  will  be  small  for  waves  that  are  well  resolved  under  the 
grids  on  both  sides  of  the  interface,  since  the  solution  discontinuity  decreases  dramatically  as  the  resolution 
of  the  scheme  increases.  This  is  further  illustrated  in  Figure  5.2  where  regions  that  satisfy  the  requirement 
on  the  resolution  (number  of  elements  per  wavelength)  so  that  the  reflection  is  2%  or  less  are  plotted  for 
a  given  grid  discontinuity  of  ratio  h^/hi.  A  value  of  7  =  0.5  is  used  in  the  calculations.  The  solid  line  is 
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Fig.  5.3.  Enlarged  numerical  phase  (a)  and  damping  (b)  errors  per  wavelength  of  propagation  for  the  physical  mode, 
is  the  numerical  wave  number  and  K  is  the  exact  wave  number,  r  =  wave  amplitude  damping  factor. 

The  dotted  lines  indicate  the  accuracy  limits  used  in  plotting  Figure  5.2. 


the  2%  reflection  boundary  for  each  given  order  of  the  scheme  as  indicated  on  the  graph.  As  we  can  see, 
when  the  ratio  h^/hi  increases,  the  requirement  on  resolution  also  increases.  It  is  interesting  to  compare 
this  requirement  with  the  resolution  requirement  placed  by  the  accuracy  criteria  of  the  scheme  had  the  grid 
been  uniformly  spaced.  Since  here  we  assume  /i2  >  hi,  the  accuracy  requirement  will  be  calculated  based  on 
/i2.  The  accuracy  boundaries  are  plotted  in  Figure  5.2  as  dotted  lines.  The  criteria  used  here  consist  of  the 
dispersion  error  27r  <  0.001  and  dissipation  error  1  —  g-‘^'^irn,(khh)/K  ^  o.OOl.  This  corresponds 

to  requiring  that  the  phase  and  damping  errors  be  less  than  10%  after  a  wave  has  been  propagated  100 
wavelengths.  Enlarged  numerical  dispersion  relations  are  plotted  in  Figure  5.3  where  the  accuracy  limits 
used  are  shown  as  dotted  lines.  Figure  5.3  indicates  that  the  uniform  grid  accuracy  constraint  is  similar 
to,  and  in  many  cases  more  stringent  than,  the  accuracy  constraint  due  to  the  abrupt  change  in  mesh  size. 
Although  the  uniform  grid  and  discontinuous  grid  error  criteria  used  here  are  somewhat  arbitrary,  we  use 
Figure  5.2  to  emphasize  the  notion  that  both  types  of  errors  follow  parallel  trends  with  respect  to  varying 
mesh  sizes  and  the  increase  of  the  resolution  of  a  scheme  leads  to  the  reduction  of  numerical  reflection  caused 
by  a  grid  discontinuity. 


20 


Fig.  6.1.  Propagation  of  a  periodic  sine  wave. 


computed  in  Table  6.1. 


The  difference  of  the  solutions  at  two  periods  shown  in  dark  lines  are 


6.  Numerical  examples.  In  this  section,  we  present  numerical  examples  that  illustrate  and  verify  the 
wave  propagation  properties  found  in  this  paper. 

6.1.  Super-accuracy  of  wave  propagation.  We  solve  the  linearized  Euler  equations  with  constant 
mean  flow  in  1-D: 


(6.1) 


du  , ,  dp 


(6.2) 


dp  , , 

dt  dx 


+ 


du 

dx 


=  0 


where  M  is  the  mean  flow  Mach  number,  u  is  the  velocity  and  p  is  the  pressure.  The  Jacobian  matrix  has 
eigenvalues  M  —  1  and  M  +  1,  which  represent  the  acoustic  wave  modes.  We  use  Legendre  polynomials  as  basis 
functions  in  our  calculation.  The  semi-discrete  equation  is  solved  by  an  optimized  4th-order  Runge-Kutta 
scheme  (LDDRK56  in  [10]). 

To  verify  the  accuracy  of  spatial  propagation,  we  consider  a  computational  domain  of  [0, 100]  and  intro¬ 
duce  an  incoming  wave 


Pin 

at  the  left  boundary  a;  =  0.  The  frequency  is  chosen  to  be  wq  =  7r/2  with  a  wavelength  Aq  =  4  in  a  mean 
flow  M  =  0.  At  the  right  boundary  x  =  100,  we  implement  the  characteristics  boundary  condition,  i.e.,  the 
exact  characteristics  flux  formula  (7  =  1)  is  used  at  the  right  boundary  of  the  last  element.  After  the  initial 
transient  has  exited  the  right  boundary,  the  computational  domain  is  filled  with  the  sine  wave.  We  then 


=  sin[a;o(a;  —  t)] 
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Table  6.1 

Solution  of  (6.1)-(6.2),  M  =  0,  using  uniform  grids.  Error  is  calculated  by  (6.4). 


7  =  1 

7  =  0.5 

p 

h 

Error  E 

order 

Error  E 

order 

1 

1.74054 

- 

1.79386 

- 

1 

0.5 

1.09166 

0.6730 

1.46813 

0.2890 

0.25 

0.197915 

2.4635 

0.344971 

2.0894 

0.125 

0.0261657 

2.9191 

0.0506057 

2.7691 

1 

0.27629 

- 

0.286715 

- 

2 

0.5 

0.010116 

4.7714 

0.00634575 

5.4976 

0.25 

0.000323692 

4.9658 

0.000172082 

5.2053 

0.125 

0.1016  X  10“^ 

4.9923 

0.5165  X  10“® 

5.0574 

1 

0.00386958 

- 

0.00381781 

- 

3 

0.5 

0.3217  X  10“^ 

6.9102 

0.4912  X  10“^ 

6.2801 

0.25 

0.2552  X  10“® 

6.9780 

0.4709  X  10“® 

6.7048 

0.125 

0.2019  X  10“® 

6.9812 

0.3964  X  10“® 

6.8919 

2 

0.0126055 

- 

0.0238191 

- 

4 

1 

0.3002  X  10“^ 

8.7137 

0.3034  X  10-^ 

9.6164 

0.5 

0.6153  X  10-^ 

8.9305 

0.3858  X  10-^ 

9.6192 

compare  the  numerical  solutions  at  the  first  period  near  a;  =  0  with  that  of  the  20th  period,  noted  by  dark 
lines  in  Figure  6.1.  Specifically,  we  measure  the  error  according  to  pressure  p  as: 


where  n  is  the  element  index  and  no  =  ^.  Table  6.1  shows  the  mesh  refinement  results  for  p  =  1  to  4. 
Since  the  local  dispersion  relation  is  accurate  to  order  2p  +  2,  the  global  error  measure  E  defined  in  (6.4) 
will  decrease  at  order  2p  +  1.  This  is  observed  in  all  the  cases. 

6.2.  Reflection  at  grid  discontinuity  and  comparison  with  eigenfunctions.  In  Figures  6.2  to 
6.3,  we  show  the  propagation  of  the  sine  wave  (6.3)  through  a  mesh  discontinuity.  Since  the  numerical  wave 
reflection  properties  are  dependent  on  the  flux  formula  used,  we  will  show  cases  with  the  exact  as  well  as 
inexact  characteristics  flux  formulas.  This  will  be  indicated  by  the  value  of  7  used  in  the  computation.  A  value 
of  I7I  =  1  indicates  exact  characteristics  splitting  while  a  value  of  I7I  ^  1  indicates  inexact  characteristics 
flux.  In  some  calculations,  a  fairly  large  grid  discontinuity  has  been  used.  This  is  to  make  reflection  errors 
more  visible  for  the  purpose  of  illustration.  In  all  the  calculations,  a  fifth-order  (p  =  4)  scheme  is  used. 

6.2.1.  Exact  characteristics  flux  formula  jyj  =  1.  In  Figures  6.2(a)  and  (b),  a  grid  discontinuity 
is  introduced  at  a;  =  50  with  the  ratio  of  grid  spacing  being  2  and  5  respectively.  The  exact  characteristics- 
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Fig.  6.2.  Propagation  of  a  periodic  sine  wave  through  a  grid  discontinuity  at  x  =  50.  M  =  0,  7  =  1  (exact  characteristics 
flux).  +  indicates  grid  points,  (a)  hi  =  1,  /i2  =  2;  (b)  hi  =  1,  /i2  =  5. 

based  flux  formula  is  used  in  this  example  with  6  =  1.  In  both  cases,  the  abrupt  change  of  element  size  causes 
no  numerical  reflection  because  the  opposite-traveling  spurious  mode  is  now  non-existent.  The  damping  of 
wave  amplitude  in  the  coarsened  grid  is  due  to  the  reduction  in  resolution  and  is  expected. 

6.2.2.  A  slow  wave  mode  I7I  =  10.  In  Figure  6.3(a),  we  show  the  solution  for  a  case  in  which  7  is 
large  (7  =  10).  This  situation  is  likely  to  occur  when  the  wave  speed  of  an  eigenmode  is  small  relative  to 
the  fastest  eigenmode  governed  by  a  given  system  of  equations.  In  Figure  6.3(a),  the  amount  of  reflection 
is  visible  since  the  grid  ratio  here  is  quite  large.  By  subtracting  out  a  calculation  with  uniform  grids  (done 
separately),  the  reflected  wave  is  extracted  and  plotted  in  Figure  6.3(b).  Inspecting  visually,  the  reflected 
wave  is  in  the  form  of  the  spurious  numerical  mode.  This  will  be  further  confirmed  when  we  compare  the 
numerical  solution  with  the  eigenfunction  formed  in  (3.20). 

To  compare  the  numerical  solution  with  the  eigenfunctions  found  in  section  4,  we  first  extract  the 
complex  coefficient  vector  from  the  numerical  solution  by  constructing 

in  each  element.  In  the  above,  T  is  the  period  of  the  sine  wave  and  to  is  an  arbitrary  time  at  which  the 
numerical  solution  has  become  time  periodic  and  v  denotes  the  solution  coefficient  vector  of  the  pressure  p. 
Then,  we  fit  this  coefficient  vector  by  a  linear  combination  of  the  eigenvectors  of  (3.18).  Specifically,  suppose 
the  eigenvectors  of  (3.18)  are  denoted  by  and  for  the  physical  and  spurious  modes,  we  try  to  And 
a  and  b  such  that 
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(6.5) 


V  =  av(P^  + 


The  coefficients  a  and  b  are  computed  by  requiring  (6.5)  be  orthogonal  to  and  .  In  other  words, 
we  “decompose”  the  numerical  solution  into  eigen-modes.  This  is  done  for  every  element  and  the  residues 
of  (6.5)  have  been  found  to  be  near  machine  zero  in  all  cases.  The  magnitudes  of  a  and  b  plotted  in  Figure 
6.3(c).  Here,  circles  indicate  the  magnitude  of  the  physical  mode,  |a|,  and  the  triangles  the  spurious  mode, 
|6|.  The  reflection  at  the  interface  at  a;  =  50  and  their  subsequent  exponential  decay  are  clearly  shown.  Also 
shown,  in  dotted  lines,  are  the  predictions  of  the  reflected  and  transmitted  waves  with  their  amplitudes  at 
the  interface  being  determined  by  (5.13)  and  (5.14).  Excellent  agreements  are  found. 

6.3.  Propagation  of  an  aconstic  pnlse  with  mean  flow.  In  the  third  example.  Figure  6.4,  we  show 
the  propagation  of  an  acoustic  pulse  in  a  mean  flow  of  Mach  number  M  =  0.8.  We  solve  (6.1)-(6.2)  using  the 
Lax-Friedrich  formula  (2.9)  with  0  =  1.  The  initial  Gaussian  profile  in  the  u  velocity  component  is  separated 
into  a  downstream  propagating  pulse,  with  speed  M  +  1,  and  an  upstream  propagating  pulse,  with  speed 
M  —  1.  Both  pulses  are  to  propagate  through  a  grid  discontinuity  of  ratio  h^/hi  =5  located  at  a;  =  30 
and  a;  =  —  30  respectively.  The  difference  in  wave  propagation  speed  results  in  two  different  upwind  factors 
7  for  the  two  pulses,  namely,  7  =  1  for  the  downstream  propagating  pulse  and  7  =  — 9  for  the  upstream 
propagating  pulse  according  to  (3.15)  and  (3.16).  For  the  right-traveling  pulse,  since  the  flux  formula  is  the 
exact  characteristics  splitting,  no  reflection  occurs  as  the  pulse  propagates  through  the  grid  discontinuity. 
For  the  left-traveling  pulse,  small  reflected  waves  are  detected  due  to  the  inexact  characteristics  flux  formula 
for  that  wave  speed.  We  note  that  the  reflected  waves  are  in  the  form  of  spurious  waves  and  decay  rapidly. 
We  emphasize  that  the  use  of  a  relatively  large  abrupt  increase  in  grid  size  is  to  make  the  reflections  more 
visible.  Indeed,  a  calculation  using  a  grid  ratio  of  2  produced  much  smaller  reflected  waves. 

7.  Conclusions.  We  have  carried  out  a  detailed  study  of  spatially  propagating  waves  in  a  discontinuous 
Galerkin  scheme  applied  to  a  system  of  linear  hyperbolic  equations.  An  eigenvalue  problem  for  the  spatially 
propagating  waves  is  formulated.  In  one  dimensional  space,  the  eigenvalue  problem  reduces  to  a  quadratic 
equation  and,  consequently,  yields  at  most  two  numerical  wave  modes  for  each  physical  wave  mode  of  the 
partial  differential  equations.  One  is  physically  significant  with  the  dispersion  error  that  decays  like 
and  the  dissipation  error  that  decays  like  locally.  The  other  numerical  mode  is  spurious.  The  spurious 

mode  becomes  non-existent  when  the  exact  characteristics  splitting  flux  formula  is  used.  Furthermore, 
reflection  and  transmission  coefficients  of  an  incident  wave  at  an  interface  of  grid  discontinuity  are  derived. 
It  is  shown  that  numerical  reflection  error  consists  of  only  the  spurious  mode  and  its  magnitude  depends 
on  the  spatial  resolution  of  the  grids  on  both  sides  of  the  interface.  Theoretical  predictions  are  verified 
with  numerical  examples.  These  predictions  should  benefit  the  design  and  application  of  the  DGM  scheme 
with  non-uniform  grids.  In  a  forthcoming  paper,  we  will  examine  the  effects  of  grid  discontinuity  in  two- 
dimensional  space. 
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Fig.  6.3.  (a)  Propagation  of  a  sine  wave  through  a  grid  diseontinuity,  =  5,  7  =  10.  (b)  Refleeted  wave,  (e) 

Decomposition  of  numerieal  solution  into  physieal  and  spurious  modes,  eireles,  physieal  mode;  triangles,  spurious  mode; 
dashed  lines,  theoretieal  predietions  of  (5.13)-(5.14). 
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Appendix. 


Al.  Kronecker  product.  Let  A  = 

{  }/  X  i; 

and  B 

—  {/^y  }nxm 

CCllB 

q:i2B 

...  O^B 

(7.1)  A®B  = 

q:2iB 

0:220 

...  02i;B 

cc/iB 

0/2  B 

...  0/A;B 

The  Kronecker  product  is  defined  as 


Inxkm 


It  is  easy  to  verify  by  direct  calculation  that  for  any  matrices  A,  B  and  vectors  x,  y,  we  have 


(7.2) 


(A  8)  B)(x  8)  y)  =  (Ax)  8)  (By) 


A2.  Polynomials  G{x)  and  H{x).  Polynomials  G{x)  and  H{x)  appeared  in  (4.2): 


G{x) 


H{x) 


p  =  l 


1  +  jX  +  ^X^ 


p  =  2 


l-lx  +  ^x^ 


1  -b  —X  -b  —x‘^  -b  —X^ 

±  -I-  gX  -I-  2qX  ^  QffX 


P  -- 

JL 

P  ■- 


--  3 
-4 
--  5 


1 


la: 


+  14 


-^x- 

200'*' 


1  ^.3  ,  1  ^4 

126"^  ^  3024'*' 

J_t3  4.  1  .p4 
99'*'  ^  1584'*' 


1 


.  4  ,  J_  2 

qX  -r  ^24- 


1  -I-  —7*  -I-  -I-  ^ -I-  ^ 

X  -r  9^.  -r  30^.  t  252^  ^  3024 

~x  H“  —x^  H“  —x^  H“ 

1  u.  -r  22^  ^99^  ^  528^ 


1  - 


-^X  +  —x^ 
11^  -r  11^ 


_J_^5 

55440^ 


It  is  straightforward  to  verify  that  H{x)/G{x)  is  exactly  the  Pade  approximation  of  e*  to  order  2p-b  2. 
This  has  been  confirmed  up  to  p  =  16  and  is  conjectured  to  be  true  for  all  p. 

A3.  Accuracy  of  eigenvalues.  In  this  appendix,  we  give  a  derivation  of  (4.5)  and  (4.6).  We  need 
only  to  consider  the  cases  when  |7|  ^  1.  The  case  for  |7|  =  1  follows  trivially  from  (4.2). 

For  convenience  of  discussion,  define 

/3  =  (_l)P+i^. 

1-7 

Then,  equation  (4.3)  becomes 


(7.3) 


A"  - 


H{iK) 


[G{iK)  ^  G{iK)  \ 

Let  the  two  roots  of  (7.3)  be  denoted  by  and  A^^^  with  their  values  at  A  =  0  as  follows, 

at  A  =  0,  A^^’)  =  1  and  A^")  =  (i. 

Now  consider  an  auxiliary  quadratic  equation: 


(7.4) 


<7^- 


+  I3e 


-rnGj-iK) 

G{iK) 


Mi-iK)  ^ 
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The  two  roots  of  (7.4)  are  easily  found  to  be 


(j(p)  =  giK 


and 


d-)  = 


G{-iK) 
G(iK)  ■ 


By  subtracting  (7.4)  from  (7.3),  we  get 


(7.5) 


A"  -cr" 


HjiK)  H{-iK) 
G{iK)  ^  G{iK) 


G{-iK) 

G(iK) 


(T  =  0. 


By  the  results  of  Appendix  A2,  we  have 


H{iK) 

G(iK) 


=  +  R{iK) 


and 


H{-iK)  ^  H{-iK)  G{-iK)  ^ 

G{iK)  G{-iK)  G{iK)  ^  ^  G{iK) 


_-iKG{-iK) 

G(iK) 


+  R{-iK) 


G{-iK) 

G(iK) 


where  R{x)  is  an  0(a;^^+^)  function  with  real  coefficients.  Then  equation  (7.5)  can  be  written  as 


(7.6) 


A" 


G{-iKy 
G{iK)  _ 


(A-cr) 


R{iK)  +  PR{-iK) 


G{-iKy 
G{iK)  _ 


A. 


For  simplicity,  let’s  define 


A(iK)  =  +  I3e 


G{iK)  ' 


B{iK)  =  R{iK)  +  l3R{-iK) 


G{-iK) 

G{iK) 


and  rewrite  (7.6)  as 


(7.7)  {X-(T)[X  +  (T-A{iK)]=B{iK)X. 

It  is  straightforward  to  verify  that,  as  K  0,  we  have 


(7.8)  A{iK)  =  [l  +  p]+^^i  (iK)  +  /i2  (iKf  +  •  •  • 

(7.9)  B{iK)  =  [1  +  /3]  i/i  {iKyP+‘^  +  +  •  •  • 

and 

(7.10)  A^^’)  +  -  A{iK)  =  [1  -  /3]  +  ii\{iK)  +  +  ■■■ 
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(7.11)  -  A{iK)  =  [p-l]+  ti'liiK)  +  n'liiKf  +  •  •  • 

where  all  the  coefficients  on  the  right  hand  sides  are  real  and  the  dots  represent  higher  order  terms  in  {iK) 
with  real  coefficients. 

Therefore, 

(i)  if  /3  ^  ±1,  it  follows  easily  from  (7.7)  that 

(7.12)  \{P)  -a(p)  =  +  ■■■ 

[-*■  P\ 

and 

(7.13)  +  . . . 

This  immediately  leads  to  (4.5)  and  (4.6). 

(ii)  if  /3  =  1,  then,  instead  of  (7.10)  and  (7.11),  we  have 

X(p)  +  (j(p)  _  A{iK)  =  ii\{iK)  +  ii'SKf  +  ■■■ 

\{s)  +  „(s)  _  A(iK)  =  ii'liiK)  +  ii'iiiKf  +  ■■■ 

which  gives 

\{P)  _(j{p)  =  2  +  ... 

and 

A(")  -  =  l-f3p^{iKYP+'^  +■■■ 

This  is  one  order  lower  than  that  of  case  (i). 

(hi)  if  /3  =  —1,  then,  instead  of  (7.9),  we  have 

B(iK)  =  iy2(iKfP+^  +  ■  ■  ■ 


and  it  follows  from  (7.7)  that 

A(p)  -cr(P)  =  l^^(iK)‘^P+3  +  . . . 

and 

AW  -crW  =  +  ••• 

This  is  one  order  higher  than  that  of  case  (ii). 

We  note,  however,  that  (ii)  or  (hi)  is  possible  only  if  7  =  0  or  7  =  oo.  Both  are  not  very  practical 
situations. 
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